Study area

This study was carried out in the National Park Berchtesgaden and its surroundings, located in the southeasternmost corner of Germany as part of the Berchtesgaden Alps, which belong to the Northern Limestone Alps. The national park is characterized by steep valleys and high mountains ranging from 600 up to 2700 m above sea level (MASL). The mean annual temperature (MAT) averaged over the last 10 years (2011–2020) varies between −0.8 and 9.7°C, depending on the elevation, with an annual precipitation of 1240–2881 mm (extracted from 1 × 1-km grid data of yearly means from the German Meteorological Service [DWD]). The higher elevations above the timberline host mainly alpine meadows, whereas coniferous forests interspersed with pastures dominate lower elevations. Alpine pastoral use has a centuries-long tradition in this region. However, some meadows are abandoned, whereas others are still traditionally managed.

In 2009, 34 grassland sites were selected, covering an elevational range from 600 up to approximately 2000 MASL (Leingärtner et al., 2014), which was the highest elevation at which a closed grass layer can be found. The sites were arranged along five transects with approximately 250 m in altitude in between. Twenty-nine sites were located inside the national park and five in the close surroundings to extend the elevational gradient towards the lowlands (Leingärtner et al., 2014). All grasslands were unused or extensively managed by extensive grazing or one late summer cut without any fertilization. In 2019, 33 of these 34 grassland sites could be revisited (Figure 1c). Owing to the national park status of most sites, management remained stable within the considered decade.

Data collection

Statistical analyses

We performed all statistical analyses using R version 4.1.1 (R Core Team, 2021). Daily temperatures were averaged per year (mean annual temperature [MAT]) from 2009 to 2019 and for the respective summer months (May–September) (mean summer temperature [MST]). We chose MST as a predictor variable because it had a better fit in our models than MAT, reflecting the sampling period and the main flight time of adult butterflies. In addition to the absolute temperature difference between the two sampling years (2009/2019), we estimated the gradual temperature change over the decade (2009–2019) by extracting the slope of linear models (MAT/MST ~ year) (Appendix S1: Figure S2), showing a steeper trend for the temperature increase over the decade than the absolute difference between the sampling years. The changes in temperature differences along the elevational gradient were analyzed using GAMs with the package mgcv version 1.8–38 (Wood, 2017), fitted with a smoothing function of elevation and setting the basis dimension k of the smoothing parameter to four (Figure 1b, Table 1). By setting k = 4, we limit the number of turns the curve is allowed to make when fitted to the data, equaling a polynomial function of maximum fourth order to avoid overfitting and keep the number of turns in an explainable dimension. We used k = 4 for all GAMs and generalized additive mixed models (GAMMs) described in the following paragraphs. We determined the local temperature lapse rate along the elevational gradient by extracting the slope of a linear model (temperature and elevation), based on mean daily temperatures from nearby climate stations, averaged over the considered decade per station (Appendix S1: Figure S3). The lapse rate and the change of temperature between or across years were used to calculate expected species range shifts and to compare them with observed shifts.

We pooled butterfly species records per site for each year and averaged the estimated flower cover per site over all visits of the year. The scale values of the vegetation survey for each host plant species were transformed into an averaged percentage cover value of all 10 subplots per site, using the Domin 2.6 transformation from Currall (1987). Species richness of butterflies and host plants was calculated with the vegan package version 2.5–7 (Oksanen et al., 2020) and analyzed using GAMMs with the mgcv package version 1.8–38 (Wood, 2017). The models were fitted with a smoothing function of elevation, year as a factor (2009, 2019), and the interaction of elevation and year, including study site as a random intercept effect and setting data family to Poisson. We used the observed species richness of butterflies for calculations because it was highly correlated (r = 0.98) with the abundance-based coverage estimator (Chao & Lee, 1992) (ACE, calculated with the package fossil V. 0.4.0 [Vavrek, 2011]) showing a high species richness saturation (73%–100%). Thus, the data seemed to reflect well the local species community. Because of the many “rare” species with few individuals, especially at high elevations with harsh weather conditions, ACE was most appropriate for validation in comparison to other indices like Chao1. It recognizes those “rare” species as such (Chao & Shen, 2004) instead of overestimating species richness at plots with many low numbers of individuals.

We used piecewise structural equation models (SEMs), provided by the piecewiseSEM package version 2.1.2 (Lefcheck, 2016), to identify (a) drivers of butterfly species richness and (b) drivers of the change in richness between years. In comparison to traditional SEMs, these models allow a wider variety of model structures and can account for various data distributions (e.g., Poisson) or relatively small sample sizes as typical issues in ecological research. We also considered the latter in model selection using the Akaike information criterion with a correction for small sample sizes (AIC_c) for decision-making. We selected four potential drivers of butterfly richness: Mean summer temperature (MST), management, adult food resources (flower cover), and number of larval food plants (host plant richness). We chose MST because the model had a higher explanatory value and better fit (lower AIC_c) than with MAT. Management per site was classified either as extensively used (1), including one late summer cut and extensive grazing, or unused (0). In addition to host plant richness, host plant cover was included initially but disregarded in the end because it added no explanatory value to the model. For the second SEM, explaining the change in butterfly richness, we included the absolute differences of all explaining variables in the SEM (e.g., ∆MST_2019–2009). Thus, sample size reduced to 32 in this latter model, and study site was no longer included as a random term. Change in management was not included because there were no significant differences in land use between years. We considered adding management itself as an explanatory variable as used in the first path model instead, which led to one alternative model to that presented in the results, where management replaced change in MST as the explanatory variable. Models were statistically not distinguishable (∆AIC_c ≤ 1.5), which was probably driven by the moderate correlation of change in MST and management (corr = −0.51, p = 0.003) in our study region. We further checked for the possible impact of the overall temperature trend over the decade (mean annual temperature 2009–2019) on species richness change as an alternative explanatory temperature variable to change in MST. Models were also not statistically different (∆AIC_c ≤ 1.5), but the explanatory value of the model got worse (R² from 0.27 to 0.21), so that we concluded that change in MST appears to be the more suitable variable to explain the change in butterfly species richness. Therefore, we here present the (best) model, including change in MST, but discuss its uncertainty.

To identify shifts of communities and species, we included butterfly species with a total count of at least 10 individuals in each year and host plant species that occurred in at least 10 subplots per year. The exact magnitude of overall butterfly shifts as well as the number of included species varied according to the selected minimum number of individuals that a species needed to be included in the analysis, ranging from 17 included species (minimum 30 individuals) to 43 included species (minimum five individuals). For host plants, it varied between 44 included species (minimum 30 occupied subplots) and 85 included species (minimum five occupied subplots). The chosen minimum numbers of individuals/subplots reflected the overall change without excluding too many species (Appendix S1: Figures S4 and S5). Hence, we included 37 out of 71 butterfly species and 70 out of 128 host plant species.

We calculated butterfly community shifts by comparing the “community elevation score” between years, following Chen et al. (2009). To calculate the score, we assigned an average occurrence elevation to each species, using the abundance-weighted mean elevation of each species in 2009. Subsequently, we averaged the elevation values of all recorded species, weighted by abundance, to obtain the final score per site for each year. We analyzed the differences in community elevation score between years using a GAMM with the mgcv package version 1.8–38 (Wood, 2017). We fitted it with a smoothing function of site elevation, year as a factor (2009, 2019), and the interaction of site elevation and year, including study site as a random intercept effect.

To evaluate mean species range shifts, we used a linear mixed effect model within the lme4 package version 1.1–27.1 (Bates et al., 2015), with the abundance-weighted mean elevation per butterfly species as response variable and year as explanatory variable, including species as a random term. The same method was applied to check mean host plant species range shifts, replacing the response variable with the mean elevation of occupied subplots per host plant species. The difference in shifts between alpine and nonalpine species was calculated using a linear model with the shifts per species explained by the mean species elevation in 2009 and the factor alpine/nonalpine. A linear model with the elevations of all recorded individuals (butterflies)/occupied subplots (host plants) as response variable and year (continuous), species (factor), and their interaction as explanatory variables was used to test for shifts of single species (Figure 5). Exact shifts, corresponding standard errors and significance of shifts were derived post hoc by calculating estimated marginal means with the emmeans package version 1.7.2 (Lenth, 2022).

Interannual changes in temperature, butterfly, and host plant richness

In 2019, MAT on our study sites was on average 0.94°C warmer than in 2009. When comparing MST, it was 0.22°C warmer in 2019 than in 2009. In both cases, the difference in temperature between the two years changed significantly along the elevational gradient with a more pronounced temperature increase at higher elevations (GAMs, p ≤ 0.001, Figure 1 [dark lines] and Table 1). Interestingly, MST at elevations up to 1000 MASL was even cooler in 2019 than in 2009. The temperature differences between 2009 and 2019 also reflected the mean gradual temperature increase, which was even higher with 1.45°C in MAT and 1.23°C in MST over this decade. Equally, the gradual temperature increase over time changed significantly along the elevational gradient, with a stronger increase at high elevations (GAMs, p ≤ 0.001, Figure 1 [light lines] and Table 1).

Butterfly abundance did not change significantly along the elevational gradient but differed significantly between years (GAMM, elevation: year p ≤ 0.001, Table 1), being significantly lower in 2019 (GAMM, intercept between years p ≤ 0.001). The distribution of species richness along the elevational gradient changed from a linear decline in 2009 to a midelevation peak at ~1200 MASL in 2019 (GAMM, elevation: year p ≤ 0.05, Figure 2a and Table 1). Only in the lowlands (<1000 MASL) did species richness decline significantly when the elevational gradient was divided into three sections (lowland [644–999 m], midelevations [1000–1499 m], highlands [1500–2034 m]) and yearly species numbers were tested against each other within sections (t-tests, lowlands: t = 2.8776, df = 17.855, p = 0.0101, midelevations: t = −0.7474, df = 21.501, p = 0.4629, highlands: t = −0.6734, df = 18.863, p = 0.5089).

Altogether, we recorded 128 potential host plant species for the complete recorded butterfly species pool, of which we found 119 in 2009 and 122 in 2019. Host plant richness peaked at low and midelevations and declined above ~1250 MASL (Figure 2b). There was no significant difference in host plant richness between 2009 and 2019 along the elevational gradient (GAMM, difference between years p > 0.1, Table 1).

Drivers of butterfly richness

Path analysis indicated that MST, management, and host plant richness had a positive effect on butterfly richness (Figure 3a,b), with MST having the greatest overall effect. MST and host plant richness affected butterfly species richness directly, whereas extensive grassland management showed an indirect positive effect on butterfly richness by increasing host plant richness. Direct, negative effects of management on butterfly richness and positive effects of flower cover were not significant but included in the best-fitting path model (Fisher's C = 1.354, df = 2, p = 0.508, N = 64) that explained 63% of the variance in butterfly species richness. There was one alternative path model (∆AIC_c ≤ 1.5) where the nonsignificant relationship between MST and host plant diversity was discarded (dashed line). In this model, only the standard estimate of the relationship from management to host plant diversity changed from 0.47** (best model) to 0.62***.

To explicitly check the mechanisms that explain the difference (∆) in species richness between 2009 and 2019, we ran a second path model (Figure 3c,d). We only included the differences between years in all explaining and response variables per site. The best-fitting path model (Fisher's C = 0.194, df = 2, p = 0.907, N = 32) including these delta values explained 27% of the change in butterfly species richness between 2009 and 2019. In this model, the change in MST was the only variable significantly affecting the change in butterfly species richness between years. There was one alternative model (∆AIC_c ≤ 1.5) where the nonsignificant relationship between ∆MST and ∆ flower cover was discarded (single arrow model: ∆MST → ∆ butterfly species richness, standard estimate 0.43, p ≤ 0.05). However, it explained a lower proportion of the variance of the ∆ in species richness (R² = 0.19). The same applied when replacing ∆MST with management or overall temperature trend over the decade (mean annual temperature 2009–2019): The path models were not statistically different from the one with ∆MST (∆AIC_c ≤ 1.5), but they explained less of the change in species richness (R² = 0.23 or R² = 0.21).

Upslope shifts of communities and species

The butterfly community elevation score, averaging the mean occurrence elevation of all detected butterfly species per site, showed a general upwards shift of species communities, as the score in 2019 was significantly lower than in 2009 (GAMM, elevation: year p ≤ 0.001, Figure 4 and Table 1). Averaging the difference in community elevation scores between years of each plot, communities shifted on average 51 ± 14 m (mean ± SD) upwards.

On average, butterfly species moved 79 ± 19 m upwards based on a comparison of the difference in mean elevation between the two years (linear mixed effect model, df = 36, p ≤ 0.001). Shifting distances decreased significantly with increasing altitude (p ≤ 0.05), but in contrast, mountainous species (lower elevational limit ≥500 MASL according to Paolucci [2013]) shifted 139 ± 57 m more than other species (linear model, df = 34, p ≤ 0.05, Appendix S1: Figure S6). At the individual species level, 14 out of 37 selected butterfly species shifted significantly along the elevational gradient between 2009 and 2019 (estimated marginal means, post hoc to linear model [df = 6312, p ≤ 0.001]), 11 upwards and three downwards (Figure 5 left). Although most mountainous species shifted (eight out of 11), only about a quarter (six out of 26) of the species occurring in lowlands changed significantly. Host plants shifted with an average of 7 ± 8 m considerably less (linear mixed effect model, df = 69, p > 0.1) and also the proportion of shifting species was lower (Figure 5 right): Only six out of 70 host plants shifted significantly upwards or downwards (estimated marginal means, post hoc to linear model [df = 8079, p ≤ 0.001]), and most of them were grasses and belonged to either the Poaceae or Cyperaceae family.

Considering the average increase of ~1.45°C in MAT over the last decade (Appendix S1: Figure S2) and a mean decrease in temperature of 0.52°C per 100 m of elevation (Appendix S1: Figure S3), the observed upslope shift of butterflies lagged well behind the expected shift of species required to track their temperature niche along the elevational gradient (~281 m). However, the reported gradual temperature increase over time also needs to be handled with care because the slope of the underlying linear model strongly depends on the years included in the model (e.g., MAT before 2009); here, we only considered a short time frame (2009–2019). Consideration of the absolute difference in MAT between 2009 and 2019 revealed that the expected shift was lower, about 179 m upwards, in accordance with an increase of ~0.93°C, which was closer to the observed shift rates.

Shifts in butterfly but not host plant richness

Although butterfly abundances were significantly lower in 2019 across the elevational gradient, the richness and composition of the regional species pool was very similar in both years. This contrasts with results from Spain, where significant overall declines in butterfly richness along elevational gradients over three decades were reported (Wilson et al., 2007) and signal the efficiency of the national park protection status in Berchtesgaden, in accordance with results from protected areas in the Italian Alps, where species richness even increased (Cerrato et al., 2019). Interestingly, we detected a clear upward shift of the maximum in butterfly species richness along the elevational gradient over time. Although in 2009 butterfly richness was highest in the lowlands, it peaked at midelevations in 2019. In our study, the hump-shaped pattern was mainly caused by the severe decline in species richness in the lowlands, whereby the majority of species were not lost but shifted to higher elevations, as generally expected under the gradually increasing temperatures detected in the region over the decade (Cerrato et al., 2019; Vitasse et al., 2021). Interestingly, within this period, the butterfly species decline in the lowlands was not compensated by species immigrating from even lower elevations or further south, similar to observations from Spain (Wilson et al., 2007). Such lowland biotic attrition has been described for the tropics, where no warmer habitats exist (Colwell et al., 2008). In temperate systems, it might indicate an unsuitable habitat matrix surrounding the national park, preventing butterfly species from following their climatic niches (Lenoir et al., 2020). Further, an intensification of anthropogenic influence might contribute to the lower species richness at lower elevations, as observed in several mountain ranges around the world (Forister et al., 2010; Gallou et al., 2017; Molina-Martínez et al., 2016). Our lowest study sites are located outside or at the outer margins of the national park. Thus, they might be more affected by land-use intensification (Viterbi et al., 2013) and, for example, higher amounts of insecticide use in the surrounding area (Brühl et al., 2021), even if local land use on the study sites stayed consistent.

In contrast to the shift in butterfly species richness patterns, we detected stable host plant richness along the elevational gradient. The hump-backed relationship between host plant species richness and elevation was also observed elsewhere (Gutiérrez et al., 2016) and matched the typical richness pattern observed for plants along elevational gradients in the Alps (Dainese & Poldini, 2012; Fontana et al., 2020; Hoiss et al., 2012). That there was no significant upward shift in richness might be due to the relatively time-limited study period. Plants react more slowly to temperature changes than mobile terrestrial insects (Vitasse et al., 2021). Over the entire last century, significant upward shifts in plant richness are visible (Lenoir et al., 2008), especially for mountainous and grassy species, which represent the majority of host plant species in our study. Also, an increase in plant richness at high elevations was observed, accelerating over time (Steinbauer et al., 2018). Thus, potential shifts might now become visible much faster than previously, which is supported by the slight tendency toward an upward shift of host plant richness already visible in our study over two years.

Temperature as main driver of butterfly richness

Butterfly species richness was directly driven by temperature (MST = mean summer temperature) and diversity of larval food resources (host plant richness), whereas host plant richness was driven by management, thereby indirectly also affecting butterfly richness. Adult food resource availability (flower cover) was not a significant predictor, but it did improve the path model, indicating its importance for butterfly species richness.

The strong positive correlation of butterfly richness with temperature is well known (Habel et al., 2021; Pellissier et al., 2013; Viterbi et al., 2013) as is the positive correlation of butterfly richness with (host) plant richness (Gutiérrez et al., 2016; Koch et al., 2013; Pellissier et al., 2013). However, in most cases, the two effects are examined separately without setting them in a common context. To analyze the drivers of butterfly richness in Britain, they were combined, finding a bigger impact of temperature in the form of growing degree-days than of host plant richness (Menéndez et al., 2007), corroborating our results. Nevertheless, land management was also not included in this study, which we found to be the main direct driver of host plant richness and, thus, indirectly driving butterfly richness. That a direct positive effect of extensive management on butterfly richness was found previously (Bubová et al., 2015; Jerrentrup et al., 2014) might be because (host) plant richness was not included as a possible intermediate factor. In addition, other factors not included in our analysis probably also exert an influence on butterfly richness, causing the remaining unexplained variance in our data. For example, structural diversity can positively affect butterfly richness (Viterbi et al., 2013), which we did not consider because data were not available for the first study year. Further, average annual precipitation has a positive effect on butterfly richness (Habel et al., 2021) but was not assessed at our study sites and is hard to predict at the site level for our region due to the complex topography. The latter might also be responsible for some temperature-correlated variance left unexplained, since microclimatic influences like topographic shading were not included, when modeling temperature for all sites, eventually causing some uncertainty in predictions.

According to the ascertained positive correlation between MST and butterfly richness, we would consequently expect an increase in butterfly richness with higher temperatures over the decade, and indeed, we identified the change in MST between years as the driver of change in butterfly richness between years as having the best explanatory power. Host plant richness did not explain the change; it also did not change between years. This finding is in accordance with other results explaining butterfly range shifts with temperature, but not with host plants (Wilson et al., 2005). However, in our study, the positive correlation between change in MST and change in butterfly richness resulted mainly from lower temperatures and similarly lower species numbers in the lowlands, whereas higher temperatures at higher elevations did not lead to an expected equally steep increase in species richness. The reason for this might be a time lag in the reaction of species richness to higher temperatures. This reaction might take more time than a decade; true abandonment of lower elevations or colonization of higher elevations would be required to observe changes in richness. This lag might also indicate that our results mirror delayed reactions to changes that had already occurred before the considered decade, suggesting a worse correlation between the observed changes in drivers and the reaction of species (richness). Further, we could not clearly separate the management and temperature effect, so that the negative response of species richness in the lowlands might also originate from management, possibly in interaction with temperature, as already shown in tropical mountains (Peters et al., 2019). Concerning the change in host plant richness, the missing link to change in MST in our path model seems to be a further indicator that the considered period might just have been too short to detect plant shifts. Thus, we speculate that butterfly richness changed to a certain degree with changing temperatures, whereas host plant richness had a bigger time lag in reacting because plants face a variety of challenges when trying to move upward, resulting in multiple types of lags (Alexander et al., 2018). Although less than in the lowlands (Bertrand et al., 2011), those lags could create an elevational gap between the peaks of butterfly and host plant richness in the future. Synthesizing other results from the Alps supports this assumption, as butterflies were the only order of studied organisms reacting in time to follow the upward-moving climatic niches (Vitasse et al., 2021).

Butterfly community and species upslope shifts

The community elevation score showed a clear upward shift of butterfly communities, as already reported from other regions in the world (Chen et al., 2011; Nieto-Sánchez et al., 2015; Wilson et al., 2007). The change was predominantly driven by losses of high-elevation species in the lowlands and range expansions of lowland species into midelevations, reflecting the observed changes in butterfly richness with losses in low elevations and an increase at midelevations. The slight inversion to higher community elevation scores in 2019 at the upper end of the elevational gradient might be a simple artifact because some alpine species monitored at the upper end of the elevational gradient might now have optimum ranges above the monitored elevational range. Alternatively, or additionally, temperature increases might also have increased the reproduction success and, thus, the abundances of alpine species. Because abundances are generally low at high altitudes, such an increase in high-altitude species abundances might outweigh the effect of upward-shifting lowland species.

Like the community elevation score, individual butterfly species moved considerably upward along the elevational gradient. The decreasing range shifts with increasing altitude are in line with the findings of previous studies (Mamantov et al., 2021; Roth et al., 2014), suggesting that species at higher altitudes might have certain traits slowing their shifts down (Cerrato et al., 2019; Mamantov et al., 2021; Rödder et al., 2021) or, conversely, that species at lower elevations are more pressured to react because they are closer to their upper thermal limits (Freeman et al., 2021). Interestingly, nevertheless, filtering for species restricted to higher elevations (only occurring at ≥500 MASL) shows that the mountainous species pool is reacting especially strongly (Figure 5 and Appendix S1: Figure S6), as also observed elsewhere (Cerrato et al., 2019; Molina-Martínez et al., 2016). It must be noted that the calculated significance and shifting distances of individual species were mainly abundance-driven, meaning that plots at lower elevations were not completely abandoned, but considerably more individuals were found in higher elevations, which is the reaction that could be expected over such a relatively short period. Some species did not show a significant response, possibly because of specific species traits. For example, widespread generalist species do not have a very specific distribution area that could shift visibly up- or downward.

Butterfly species in our study shifted stronger per decade than in other studies over 20–35 years (Forister et al., 2010; Molina-Martínez et al., 2016; Wilson et al., 2005). This might reflect the accelerating speed of climate change within the last decade. However, the shifts we found were still considerably below the predicted shifts to follow the climatic niches along the elevational gradient according to the observed temperature increase. Although we might have overestimated the predicted shifts owing to the limited considered time frame, a lag between observed and predicted shifts would also persist when assuming lower temperature increases in the Alps, as suggested by observations of the last 50 years (Vitasse et al., 2021). Lags were also reported in other studies (Chen et al., 2011; Devictor et al., 2012; Lenoir et al., 2020; Nieto-Sánchez et al., 2015) and appear more pronounced in temperate species communities, where species might partly benefit from higher temperatures (Freeman & Class Freeman, 2014). In temperate regions, species are still further from their upper thermal limits than in the tropics (Ghalambor, 2006), increasing the chances that species might just expand their ranges into higher elevations instead of moving their elevational range upward. Previous observations from the Alps strengthen this hypothesis, detecting stronger shifts in upper range limits than in the mean elevation of occurrence (Vitasse et al., 2021).

Nevertheless, the observed butterfly uphill shifts in mean elevation were still considerable (79 ± 19 m) and slightly above other recent observations from the Alps detecting ~40–50 m uphill shift/decade (Rödder et al., 2021; Roth et al., 2014; Vitasse et al., 2021). Thus, even if there might be an especially pronounced temporary response to the extraordinary hot year 2019 reflected in our data, a considerable species shift was found consistently throughout the Alps. At the same time, host plant species showed nearly no upward movement (7 ± 8 m), coinciding with the observations of the richness pattern and with other recent alpine studies showing similar range shifts for plants (Roth et al., 2014; Vitasse et al., 2021) when averaging shifts of alpine and herbaceous plants in the second study. The vast majority of host plants showed more uphill movement than grasses, which can be disregarded when looking at butterflies, because butterfly larvae feeding on grasses are mostly specialized at the family level, not on single plant species. Thus, moving away from a Poaceae or Cyperaceae species has no effect on them because they will find other related species. This fits together with our findings of a missing impact of host plant richness on changes in butterfly richness between years. Thus, when butterflies continue to follow their temperature niches uphill, this might form a gap between butterfly and host plant distributions, forcing them to adapt. A possible way for butterflies to avoid shifting away from their host plants is by adapting their phenology instead (Vitasse et al., 2021). The connection between host plant availability and this kind of adaptation has been demonstrated, since especially butterflies more dependent on certain larval host plants were observed to exhibit this strategy (Diamond et al., 2011). Otherwise, butterflies might also compensate for faster upward shifts broadening their larval diet breadth by including new plant species in their diet, as plants at higher elevations are less resistant to herbivores (Pellissier et al., 2012).